function matrix_simulation_expression()
    %[t1, t2, g, f] = createSymbols();

    [t1, t2] = createSymbols();

    S0 = 1; S1 = 0; S2 = 0; S3 = 0;
    incident_beam = [S0; S1; S2; S3];

    t2 = 0;g = 0; f = 90;

    rotated_polarizer = 0.5 * [1, cosd(2*g)*cosd(2*t1), cosd(2*g)*sind(2*t1), 0;
        cosd(2*g)*cosd(2*t1), cosd(2*t1)^2 + sind(2*g)*sind(2*t1)^2, (1 - sind(2*g))*sind(2*t1)*cosd(2*t1), 0;
        cosd(2*g)*sind(2*t1), (1 - sind(2*g))*sind(2*t1)*cosd(2*t1), sind(2*t1)^2 + sind(2*g)*cosd(2*t1)^2, 0;
        0, 0, 0, sind(2*g)];
    rotated_retarder = [1, 0, 0, 0;
        0, (cosd(2*t2))^2 + cosd(f)*sind(2*t2)^2, (1 - cosd(f))*sind(2*t2)*cosd(2*t2), -sind(f)*sind(2*t2);
        0, (1 - cosd(f))*sind(2*t2)*cosd(2*t2), sind(2*t2)^2 + cosd(f)*cosd(2*t2)^2, sind(f)*cosd(2*t2);
        0, sind(f)*sind(2*t2), -sind(f)*cosd(2*t2), cosd(f)];

    emerging_beam = simplify(rotated_retarder * rotated_polarizer * incident_beam);
    %emerging_beam = simplify(rotated_polarizer * incident_beam);

    fprintf('Emerging Beam Expressions:\n');
    fprintf('S0: %s\n', char(emerging_beam(1)));
    fprintf('S1: %s\n', char(emerging_beam(2)));
    fprintf('S2: %s\n', char(emerging_beam(3)));
    fprintf('S3: %s\n', char(emerging_beam(4)));
end


function [t1, t2] = createSymbols()
    syms t1 t2 real
end


%function [t1, t2, g, f] = createSymbols()
    %syms t1 t2 g f real
%end
